require(multtest)
ifst <- seq(from=0.001,to=0.99,length.out=100)
itheta <- 1/ifst -1
lik.theta.fit <- matrix(0,nrow=nrow(grouse.mar.df3[grouse.mar.df3$p.hat.all >0 &grouse.mar.df3$p.hat.all < 1,]),ncol=100,dimnames=list(grouse.mar.df3$id[grouse.mar.df3$p.hat.all >0 &grouse.mar.df3$p.hat.all < 1],paste("theta",itheta,sep="_")))


#grouse.mar.df3[1,]
#LikBetaBinomialPops(depth=grouse.mar.df3[1,c('pop1.depth','pop1.depth','pop1.depth')],
#	count=grouse.mar.df3[1,c('pop1.var.reads','pop1.var.reads','pop1.var.reads')],
#	p.hat=grouse.mar.df3[1,'p.hat.all'], theta= 0.7)

lik.theta.list2 <- vector(mode='list',length=3)
names(lik.theta.list2) <- c('all.pops','pop1v2','pop12v3')

for(i in 1:100){
	lik.theta.fit[,i] <- apply(grouse.mar.df3[grouse.mar.df3$p.hat.all >0 &grouse.mar.df3$p.hat.all < 1,c('pop1.var.reads','pop2.var.reads','pop3.var.reads', 'pop1.depth','pop2.depth','pop3.depth','p.hat.all')],1,function(X,theta){
		LikBetaBinomialPops(depth=X[c('pop1.depth','pop2.depth','pop3.depth')],
		count=X[c('pop1.var.reads','pop2.var.reads','pop3.var.reads')],
		p.hat= X['p.hat.all'],theta=theta)
		},theta=itheta[i])

	}

null.which <- which(colSums(lik.theta.fit,na.rm=TRUE)==max(colSums(lik.theta.fit,na.rm=TRUE))) # fst = 0.071
ifst[null.which]

#cut and paste from previous analysis: beta-binomial method
require(multtest)
BBlikTest.mar.df <- data.frame(null.lik = lik.theta.fit[,null.which],max.lik = rep(0,nrow(lik.theta.fit)))
BBlikTest.mar.df$max.lik <- apply(lik.theta.fit,1,max)
BBlikTest.mar.df$chi.sq <- 2*(BBlikTest.mar.df$max.lik-BBlikTest.mar.df$null.lik)
BBlikTest.mar.df$p.val <- 1-pchisq(BBlikTest.mar.df$chi.sq,1)

tmp <-  mt.rawp2adjp(BBlikTest.mar.df$p.val,'Bonferroni')
BBlikTest.mar.df$bonferroni <- tmp$adjp[order(tmp$index),'Bonferroni']

tmp <-  mt.rawp2adjp(BBlikTest.mar.df$p.val,'Hochberg')
BBlikTest.mar.df$hochberg <- tmp$adjp[order(tmp$index),'Hochberg']

#repeat with a dodgy method to pool snps at each locus

BBlikTest.mar.pooled <- aggregate(lik.theta.fit,by=list(grouse.mar.df3$cDNA[grouse.mar.df3$p.hat.all >0 &grouse.mar.df3$p.hat.all < 1]),sum)
BBlikTest.mar.pooled.df <- data.frame(cDNA=BBlikTest.mar.pooled[,1],null.lik = BBlikTest.mar.pooled[,null.which+1],
	max.lik = rep(0,nrow(BBlikTest.mar.pooled)))
BBlikTest.mar.pooled.df$cDNA <- as.character(BBlikTest.mar.pooled.df$cDNA)	
BBlikTest.mar.pooled.df$max.lik <- apply(BBlikTest.mar.pooled[,-1],1,max)
BBlikTest.mar.pooled.df$chi.sq <- 2*(BBlikTest.mar.pooled.df$max.lik-BBlikTest.mar.pooled.df$null.lik)
BBlikTest.mar.pooled.df$p.val <- 1-pchisq(BBlikTest.mar.pooled.df$chi.sq,1)

tmp <-  mt.rawp2adjp(BBlikTest.mar.pooled.df$p.val,'Bonferroni')
BBlikTest.mar.pooled.df$bonferroni <- tmp$adjp[order(tmp$index),'Bonferroni']
tmp <-  mt.rawp2adjp(BBlikTest.mar.pooled.df$p.val,'Hochberg')
BBlikTest.mar.pooled.df$hochberg <- tmp$adjp[order(tmp$index),'Hochberg']
BBlikTest.mar.pooled.df$fst <- 0

BBlikTest.mar.pooled.df$fst <- apply(BBlikTest.mar.pooled[,-1],1,function(X){
	fst.seq <- seq(from=0.001,to=0.99,length.out=100) #remember to chage this is needed
	fst.seq[which.max(X)]
	})
BBlikTest.mar.pooled.df$fst <- as.numeric(BBlikTest.mar.pooled.df$fst)
BBlikTest.mar.pooled.df <- BBlikTest.mar.pooled.df[!is.na(BBlikTest.mar.pooled.df$null.lik),]

#some big differences with population 3 seen, essentially private alleles

lik.theta.list2[['all.pops']] <- list(lik.theta.fit=lik.theta.fit,ifst=ifst,itheta=itheta, BBlikTest.df=BBlikTest.mar.df,BBlikTest.pooled.df=BBlikTest.mar.pooled.df,null.which=null.which, null.fst=ifst[null.which],cDNA=grouse.mar.df3$cDNA[grouse.mar.df3$p.hat.all >0 &grouse.mar.df3$p.hat.all < 1])


#try to repeat analysis just on pop1 vs pop2

ifst <- seq(from=0.001,to=0.5,length.out=100)
itheta <- 1/ifst -1
lik.theta.fit <- matrix(0,nrow=nrow(grouse.mar.df3[grouse.mar.df3$p.hat.1v2 >0 &grouse.mar.df3$p.hat.1v2 < 1,]),ncol=100,dimnames=list(grouse.mar.df3$id[grouse.mar.df3$p.hat.1v2 >0 &grouse.mar.df3$p.hat.1v2 < 1],paste("theta",itheta,sep="_")))


#grouse.mar.df3[1,]
#LikBetaBinomialPops(depth=grouse.mar.df3[1,c('pop1.depth','pop1.depth','pop1.depth')],
#	count=grouse.mar.df3[1,c('pop1.var.reads','pop1.var.reads','pop1.var.reads')],
#	p.hat=grouse.mar.df3[1,'p.hat.all'], theta= 0.7)

for(i in 1:100){
	lik.theta.fit[,i] <- apply(grouse.mar.df3[grouse.mar.df3$p.hat.1v2 >0 &grouse.mar.df3$p.hat.1v2 < 1,c('pop1.var.reads','pop2.var.reads', 'pop1.depth','pop2.depth','p.hat.1v2')],1,function(X,theta){
		LikBetaBinomialPops(depth=X[c('pop1.depth','pop2.depth')],
		count=X[c('pop1.var.reads','pop2.var.reads')],
		p.hat= X['p.hat.1v2'],theta=theta)
		},theta=itheta[i])

	}

null.which <- which(colSums(lik.theta.fit,na.rm=TRUE)==max(colSums(lik.theta.fit,na.rm=TRUE))) # fst = 
ifst[null.which] # 0.011


BBlikTest.mar.df1v2 <- data.frame(null.lik = lik.theta.fit[,null.which],max.lik = rep(0,nrow(lik.theta.fit)))
BBlikTest.mar.df1v2$max.lik <- apply(lik.theta.fit,1,max)
BBlikTest.mar.df1v2$chi.sq <- 2*(BBlikTest.mar.df1v2$max.lik-BBlikTest.mar.df1v2$null.lik)
BBlikTest.mar.df1v2$p.val <- 1-pchisq(BBlikTest.mar.df1v2$chi.sq,1)

tmp <-  mt.rawp2adjp(BBlikTest.mar.df1v2$p.val,'Bonferroni')
BBlikTest.mar.df1v2$bonferroni <- tmp$adjp[order(tmp$index),'Bonferroni']

tmp <-  mt.rawp2adjp(BBlikTest.mar.df1v2$p.val,'Hochberg')
BBlikTest.mar.df1v2$hochberg <- tmp$adjp[order(tmp$index),'Hochberg']

#repeat with a dodgy method to pool snps at each locus

BBlikTest.mar.pooled.1v2 <- aggregate(lik.theta.fit,by=list(grouse.mar.df3$cDNA[grouse.mar.df3$p.hat.1v2 >0 &grouse.mar.df3$p.hat.1v2 < 1]),sum)
BBlikTest.mar.pooled.df1v2 <- data.frame(cDNA=BBlikTest.mar.pooled.1v2[,1],null.lik = BBlikTest.mar.pooled.1v2[,null.which+1],
	max.lik = rep(0,nrow(BBlikTest.mar.pooled.1v2)))
BBlikTest.mar.pooled.df1v2$cDNA <- as.character(BBlikTest.mar.pooled.df1v2$cDNA)	
BBlikTest.mar.pooled.df1v2$max.lik <- apply(BBlikTest.mar.pooled.1v2[,-1],1,max)
BBlikTest.mar.pooled.df1v2$chi.sq <- 2*(BBlikTest.mar.pooled.df1v2$max.lik-BBlikTest.mar.pooled.df1v2$null.lik)
BBlikTest.mar.pooled.df1v2$p.val <- 1-pchisq(BBlikTest.mar.pooled.df1v2$chi.sq,1)

tmp <-  mt.rawp2adjp(BBlikTest.mar.pooled.df1v2$p.val,'Bonferroni')
BBlikTest.mar.pooled.df1v2$bonferroni <- tmp$adjp[order(tmp$index),'Bonferroni']
tmp <-  mt.rawp2adjp(BBlikTest.mar.pooled.df1v2$p.val,'Hochberg')
BBlikTest.mar.pooled.df1v2$hochberg <- tmp$adjp[order(tmp$index),'Hochberg']
BBlikTest.mar.pooled.df1v2$fst <- 0

BBlikTest.mar.pooled.df1v2$fst <- apply(BBlikTest.mar.pooled.1v2[,-1],1,function(X){
	fst.seq <- seq(from=0.001,to=0.5,length.out=100) #remember to chage this is needed
	fst.seq[which.max(X)]
	})
BBlikTest.mar.pooled.df1v2$fst <- as.numeric(BBlikTest.mar.pooled.df1v2$fst)
BBlikTest.mar.pooled.df1v2 <- BBlikTest.mar.pooled.df1v2[!is.na(BBlikTest.mar.pooled.df1v2$null.lik),]

lik.theta.list2[['pop1v2']] <- list(lik.theta.fit=lik.theta.fit,ifst=ifst,itheta=itheta, BBlikTest.df=BBlikTest.mar.df1v2,BBlikTest.pooled.df=BBlikTest.mar.pooled.df1v2,null.which=null.which, null.fst=ifst[null.which],cDNA=grouse.mar.df3$cDNA[grouse.mar.df3$p.hat.1v2 >0 &grouse.mar.df3$p.hat.1v2 < 1])


#should also contrast pop12 vs pop3

ifst <- seq(from=0.001,to=0.99,length.out=100)
itheta <- 1/ifst -1
lik.theta.fit <- matrix(0,nrow=nrow(grouse.mar.df3[grouse.mar.df3$p.hat.12v3 >0 &grouse.mar.df3$p.hat.12v3 < 1,]),ncol=100,dimnames=list(grouse.mar.df3$id[grouse.mar.df3$p.hat.12v3 >0 &grouse.mar.df3$p.hat.12v3 < 1],paste("theta",itheta,sep="_")))


for(i in 1:100){
	lik.theta.fit[,i] <- apply(grouse.mar.df3[grouse.mar.df3$p.hat.12v3 >0 &grouse.mar.df3$p.hat.12v3 < 1,c('pop12.var.reads','pop3.var.reads', 'pop12.depth','pop3.depth','p.hat.12v3')],1,function(X,theta){
		LikBetaBinomialPops(depth=X[c('pop12.depth','pop3.depth')],
		count=X[c('pop12.var.reads','pop3.var.reads')],
		p.hat= X['p.hat.12v3'],theta=theta)
		},theta=itheta[i])

	}

null.which <- which(colSums(lik.theta.fit,na.rm=TRUE)==max(colSums(lik.theta.fit,na.rm=TRUE))) # fst = 
ifst[null.which] # 0.0609394


BBlikTest.mar.df12v3 <- data.frame(null.lik = lik.theta.fit[,null.which],max.lik = rep(0,nrow(grouse.mar.df3[grouse.mar.df3$p.hat.12v3 >0 &grouse.mar.df3$p.hat.12v3 < 1,])))
BBlikTest.mar.df12v3$max.lik <- apply(lik.theta.fit,1,max)
BBlikTest.mar.df12v3$chi.sq <- 2*(BBlikTest.mar.df12v3$max.lik-BBlikTest.mar.df12v3$null.lik)
BBlikTest.mar.df12v3$p.val <- 1-pchisq(BBlikTest.mar.df12v3$chi.sq,1)

tmp <-  mt.rawp2adjp(BBlikTest.mar.df12v3$p.val,'Bonferroni')
BBlikTest.mar.df12v3$bonferroni <- tmp$adjp[order(tmp$index),'Bonferroni']

tmp <-  mt.rawp2adjp(BBlikTest.mar.df12v3$p.val,'Hochberg')
BBlikTest.mar.df12v3$hochberg <- tmp$adjp[order(tmp$index),'Hochberg']

#repeat with a dodgy method to pool snps at each locus

BBlikTest.mar.pooled.12v3 <- aggregate(lik.theta.fit,by=list(grouse.mar.df3$cDNA[grouse.mar.df3$p.hat.12v3 >0 &grouse.mar.df3$p.hat.12v3 < 1]),sum)
BBlikTest.mar.pooled.df12v3 <- data.frame(cDNA=BBlikTest.mar.pooled.12v3[,1],null.lik = BBlikTest.mar.pooled.12v3[,null.which+1],
	max.lik = rep(0,nrow(BBlikTest.mar.pooled.12v3)))
BBlikTest.mar.pooled.df12v3$cDNA <- as.character(BBlikTest.mar.pooled.df12v3$cDNA)	
BBlikTest.mar.pooled.df12v3$max.lik <- apply(BBlikTest.mar.pooled.12v3[,-1],1,max)
BBlikTest.mar.pooled.df12v3$chi.sq <- 2*(BBlikTest.mar.pooled.df12v3$max.lik-BBlikTest.mar.pooled.df12v3$null.lik)
BBlikTest.mar.pooled.df12v3$p.val <- 1-pchisq(BBlikTest.mar.pooled.df12v3$chi.sq,1)

tmp <-  mt.rawp2adjp(BBlikTest.mar.pooled.df12v3$p.val,'Bonferroni')
BBlikTest.mar.pooled.df12v3$bonferroni <- tmp$adjp[order(tmp$index),'Bonferroni']
tmp <-  mt.rawp2adjp(BBlikTest.mar.pooled.df12v3$p.val,'Hochberg')
BBlikTest.mar.pooled.df12v3$hochberg <- tmp$adjp[order(tmp$index),'Hochberg']
BBlikTest.mar.pooled.df12v3$fst <- 0

BBlikTest.mar.pooled.df12v3$fst <- apply(BBlikTest.mar.pooled.12v3[,-1],1,function(X){
	fst.seq <- seq(from=0.001,to=0.99,length.out=100) #remember to chage this is needed
	fst.seq[which.max(X)]
	})
BBlikTest.mar.pooled.df12v3$fst <- as.numeric(BBlikTest.mar.pooled.df12v3$fst)
BBlikTest.mar.pooled.df12v3 <- BBlikTest.mar.pooled.df12v3[!is.na(BBlikTest.mar.pooled.df12v3$null.lik),]


lik.theta.list2[['pop12v3']] <- list(lik.theta.fit=lik.theta.fit,ifst=ifst,itheta=itheta, BBlikTest.df=BBlikTest.mar.df12v3,BBlikTest.pooled.df=BBlikTest.mar.pooled.df12v3,null.which=null.which, null.fst=ifst[null.which],cDNA=grouse.mar.df3$cDNA[grouse.mar.df3$p.hat.12v3 >0 &grouse.mar.df3$p.hat.12v3 < 1])

